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Abstract 

In the natural world, there are many insect species whose individual 
members have a life history that takes them through two stages, im- 
mature and mature. Moreover, the rates of survival, development, and 
reproduction almost always depend on age, size, or development stage. 
Keeping this in mind, in this paper, a three species crop-pest-natural 
enemy food chain model with two stages for natural enemies is investi- 
gated. Using characteristic equations, a set of sufficient conditions for 
local asymptotic stability of all the feasible equilibria is obtained. More- 
over, using approach as in teeretta and KTiand . I2nn2l ). the possibility 



of the existence of a Hopf bifurcation for the interior equilibrium with 
respect to maturation delay is explored, which shows that the matura- 
tion delay plays an important role in the dynamical behavior of three 
species system. Also obtain some threshold values of maturation delay 
for the stability-switching of the particular system. In succession, using 
the normal form theory and center manifold argument, we derive the 
explicit formulas which determine the stability and direction of bifur- 
cating periodic solutions. Finally, a numerical simulation for supporting 
the theoretical analysis is given. 



Keywords: Food chain, maturation delay, stability- switch, Hopf bifurcation, 
chaos. 



1 Introduction 



It is a well known fact that pest is a harmful insec t and its outbreak of- 



ten cause ser i ous e cological and economic problems (IKaminska et al.l . 12004 



Weaver et al.l . Il992l ). Evidence indicates that annually the pests cause 25% 



loss in rice, 5-10% i n wheat, 30% in puls e s, 35% in oilseeds, 20% in sugarcane 
and 50% in cotton (jPhaliwal and Aroral . Il996[ ). Now a days, many pest con- 
trol meth ods are available, such as biological, cul t ural, physical and chemical 



methodsj Frang . 1 19611 : IVan Lenteren and Woetsl . Il988l : IVincent et al.l . 12003 



Hoytl . Il969l ). However farmers mostly use pesticides to control pests because 
of its efficiency and convenience. The chemical pesticide kills not only pests 
but it also kills their natural enemies. Actually, when pests are caught or poi- 
soned to a large extent, their natural enemies become extinct due to no food 
and afterwards, when the pesticide intensity decreases, then the pests increase 
rapidly. Due to this reason, chemical control has become challenged. Further- 
more, the common practice proves that long-term adopting chemical control 
may give rise to disastrous results, for example, environmental contamination, 
toxicosis of the man and animals and so on. Thus the pesticide pollution is also 
recognized as a majo r heal t h hazard to hurnan be i ngs and to natural ene mies 



flButleii . Il969l : iDahall . Il995l : iKaminska et al.l . 12004 iKishimba et al.l . l2004f ). On 



the other hand, it is a well known that the biological control method is harm- 
less to human, animal and environment. Biological control is generally used 
to control a particular pest using a chosen living organism; this chosen or- 
ganism might be a predator, parasite or disease which attacks on the harmful 
insect pest. The last few years have seen an sudden increase of int erest in the 
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Moreover, effective use of biological control often requires a good under- 
standing of the biology of the pest, its natural enemies and their interaction, 
as well as the ability to identify various life stages of relevant insects in the 
crops. Again, the formulations of stage structured population models are 
the recognition that individuals of many species have life-histories composed 
of a sequence of stages within which their characteristics are broadly simi- 
lar to those of other individuals in the same stage and totally different, from 
those of individuals in other stages. In insect population, such stages are 
particularly easy to recognize, being separated by short events such as moult 



(Barclav and van den Driessche. 
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Bence and Nisbet. 
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Gurnev et al. 
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Gurnev and Nisbet 
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1983 




Landahl and Hansen. 


1975 



Wood et al.l . Il989l ). but the real interest comes into pict ure on the stage struc- 



tured models after the work of Aiello and Freedman (lAiello and Freedmanl . 



I990I ). They proposed a single species model with stage structure assuming an 



average age to maturity ( constant time delay) which reflecting a de- 

layed birth of immature and a reduced survival of immature to their maturity. 



The model is as follows: 



dx(t) 
dt 



Py{t)-rx{t)-Pe-'^y{t-r), 



dyjt) 
dt 



'1.1] 



where x(t) and y(t) represent the immature and mature populations densities, 
respectively. Here, it is assumed that at any time t > 0, growth rate of 
immature population is proportional to the existing mature population with 
proportionality constant f3; the death rate of immature population is r; the 
death rate of mature population is proportional to the square of the population 
with the proportionality constant rj. The term Pe~'^'^y{t — r) represents the 
immature who were born at time t — r and survive at time t, therefore it 
represents the transformation of immature to mature, where r represent a 
constant time to maturity. All the parameters r, /3, r and r] are positive 
constants. 

Further, the single species model (11.11) is extended by many researchers 
i nto different kinds of stage-structured models and obtained significant results 



flAie: 



1994 



Huo et al 



2001 



lo et al.l.ll992l:ICao et al.l.ll992l: iFreedman and Wul . Il99ll : iFreedman et al. 



Magnusson . Il999h . Recently, many authors studied 



different kinds of predator-prey system with division of the predators into im- 
mat ure and mat u re class and a g ood number of research has been c arried 
out flWang et al l. I2OOII: IXu et all. [20oi: lOao et Zl. I2OO8I: iMagnussonl. Il999| : 



Sun et al.l . l2009l : lQu and Weil . boo Au and Huangl . l2010l : Tsatio and Takeuchil, 
2003[) . One interested model is suggested by Satio and Takeuchi (ISatio and Takeuchi 
2OO3I ). They considered two life stages for predator and proposed the following 
predator-prey model: 



x{t) = x{t) (ri - aiix{t) - arsvit)) 



Y{t) = -r2F(t) + a3ix{t)y{t) - a3ie-'-'^x{t - T)y{t - r), 
yit) = -r^y\t) + a3ie-^2^x(t - T)y{t - r). 



;i.2) 



where x{t) is population density of prey, Y{t) and y{t) denote the densities of 
immature and mature predator population, respectively; r represent a constant 
time to maturity for predator; ais is the per capita rate of predation; 031 is 
the conversion rate and all other parameters have the similar meaning as in 

Furthe rmore, three species food chain models are investigated by many re- 

searchers (Freedman and Waltrnanl. 119771: IFreedman and Sd.ll985tlHastings and Powell 



1991: Freedman and Ruan 



1992 



McCann and Yodzis. 1995: Boer et al.. 1999 



Li and Kuangl . l2000l:lHsu et al.l. 120031) . One of noteworthy contribution is given 
by Kuang et al. (IHsu et al.l . 120031 ). they considered a three trophic food chain 



model for plant-pest-natural enemy, but they ignored the stage structure phe- 
nomena of species. 

The aim of this paper is to study a crop-pest- natural enemy model with 
two life stages of natural enemy. The paper is organized as follows: in section 
|2l model development is discussed, the positivity and boundedness are estab- 
lished in the section [Hi In section SJ all the feasible equilibria and their local 
stability behavior are studied. The stability and direction of Hopf bifurcation 
is analyzed in section [51 Further, in section |6l a set of numerical simulations 
is given to verify all the major analytical findings. Finally, conclusions for this 
paper are given in the last section. 



2 Proposed Mathematical Model 



In this section, our main aim is to propose a mathematical model for the in- 
teraction of plant-pest-natural enemy. Since plant-hoppers are serious pests 
for the rice crops and these are suppressed by Lycosa tarantula and other spi- 
ders, it is well doc umented in a report of the Indian Council of Agricultural 
Research (ICAR) (IBirthal and Sharmal . |2004| ). The tarantula species has two 
major life stages, namely, immature and mature; only mature population can 
harvest the pest and reproduce a new offspring. In modelling process, we 
assume that x{t), y{t), zi{t) and Z2{t) are densities of crop, pest, immature 
and mature natural enemy at time t, respectively. The parameters ai, bi are 
respectively the intrinsic growth rate and overcrowding rate of crops; Ci, ai 
are per capita predation rate of crop by the pest and the corresponding growth 
rate of pest, respectively. The parameters C2 is per capita predation rate of 
pest by the natural enemy and ^2 is the corresponding growth rate of mature 
natural enemy. Here, di, ^2, d^ are the natural death rates of pests, immature 
and mature natural enemies, respectively. Further, r is the maturation delay 
from immature to mature natural enemies, the term a2e~'^''^y{t — T)z2{t — t) 
represents the transformation of immature to mature population. Keeping 
this biological situation in rnind and motivated from the modelling ideas of 
( lAiello and Freedmanl . Il990l : ISatio and Takeuchil . |2003[ ). in this paper, we pro- 
pose a three species stage structured crop-pest-natural enemy model as follows: 



dx(t) 
dt ' 

dyit) 
dt 

dziit) 

dt 
dz2{t) 

dt 



x{t) (ai - hix{t) - ciy{t)) , 
y{t) {aix{t) -di- C2Z2it)) , 

= a2y{t)z2{t) - d2Zi{t) - a2e~'^^''y{t - T)z2{t - r), 
= a2e-'^^^y{t - T)z2{t - r) - d^Z2{t), 



(2.1) 



The model completes with the following set of initial conditions: 



Z2ie) = MO), MO) > 0, ^.(0) >o, ee [-r,o], z = 1,2, (2.2) 

where (0i, 02, V'l, ^^2) G C ([— r, 0], i?^) , the Banach space of continuous func- 
tions mapping on the interval [— t, 0] into R^. For continuity of the initial 
conditions, we further require 

^0 

MQ) = J a2Ms)Ms)e''''ds, (2.3) 

where ■01 (0) represents the accumulated survivors of those natural enemy mem- 
bers who were born between — r and 0. In the next section, we will discuss 
the positivity and boundedness of solutions of the system (12.11) with initial 
conditions (O and (1231) . 

3 Positivity and boundedness 

Positivity means that the species is persistent and boundedness implies a nat- 
ural restriction. Therefore, we can state and prove the following lemmas for 
the positivity and boundedness of solutions of the system (12.11) : 

Lemma 3.1. The solutions of system Ii2.1\) with initial conditions ^2.2^ and 
Ii2.3\) are positive, for all t > 0. 

Proof. Let {x{t),y{t), zi{t), Z2{t)) be a solution of system (12. ip with initial 
conditions (12. 2p and (12.30 . Let us first consider Z2{t) for t G [0, r]. Noting that 
4>2{G) > 0, '02(6') > for 6 G [— r, 0], we obtain from the fourth equation of 
system (12.10 that 

= a2e-'^^"02(t - r)02(t - r) - c/3^2(t) > -4^2(t). 
It thus follows for t G [0, r], 
Z2{t) > Z2{0)e-'^^' > 0. 
For t G [0, r], the second equation of the system (12.11) can be rewritten as 

^ > y{t) (-rfi - cMt)) . 

A standard comparison argument shows that for t G [0, r]. 



y{t) > y{0) exp (^^ {-di - ciZ2is)ds)'^ > 0. 



Similarly, it follows from the first equation of system f l2.ip that for t G [0, r], 

> -bix\t) - c,x{t)y{t), 

which evidences that 

x(0)exp (-ci J^y{s)ds) 

x{t) > -j^ — ^ > 0. 

1 + 6ix(0) Jq exp (-ci /g y{u)du) 

In a similar way, we can treat the intervals [r, 2r], . . . , [nr, (n + l)r], G A^. 
Thus by induction, we establish that x{t) > 0, y{t) > and Z2{t) > for all 
t > 0. 

Finally from third equation of (12.11) and fl2.3p . we have 
zi{t)= [ a2y{s)z2{s)e-'''^'-'Us. 

Jt-T 

Therefore the positivity of Zi{t) follows. □ 

Lemma 3.2. The solutions of the system /12. 1\) with initial conditions ^2.2^) 
and Ii2.3\) are bounded. 

Proof. Let V{t) = aia2x{t) + Cia2y(t) + ciC2Zi(t) + ciC2Z2{t), calculating the 
derivative of V{t) with respect to t along the positive solution of the system 
(12.11) . we have 

V{t) = aiaia2x{t) — biaia2X^ (t) — Cidia2y{t) — ciC2d2Zi{t) — ciC3dsZ2{t) . 
Taking p = min{ai, di, d2, d^}, we obtain that 

V{t) +pV{t) < 2aiaia2x{t) — biaia2x'^ (t) . 
Hence there exists a positive constant K = a{aia2/hi such that 

V{t)+pV{t) <K, 
thus, we get 

V{t) < (V(0) - K/p) e-P*. 

Therefore, V(t) is bounded, i.e., each solution of the system (12. ip is bounded. 

□ 

In the next section, we will investigate the feasible equilibrium of the system 
(12. ip and study their stability. 



4 Nonnegative equilibria and their local sta- 
bility 

In this section, our main objective is to investigate the local behavior of all 
feasible equilibria and existence of a Hopf bifurcation at interior equilibrium. 
The equation for the variable Zi in the system (12. ip can be rewritten as 

= a2y{t)z2{t)-d2Zi{t)-a2e-'''^y{t-T)z2{t-T) 
:= -d2Zi{t) + / {y{t), Z2{t), y{t - r), Z2{t - r)) , 

if y{t), Z2{t) are bounded and y{t) — y* , Z2 Z2 as t — )■ 00, then zi{t) — )■ 
/ (y*, y*, Z2) /(i2 as t — )■ 00, i.e., the asymptotic behavior of Zi is completely 
dependent on y{t) and Z2{t). Hence, the asymptotic behavior of our proposed 
model will remain the same with the following reduced system: 

"^^^^^ - x{t) (ai - bix{t) - ciy{t)) , 



dt 

dyit) 
dt 



y{t) {aix{t) - c/i - C2Z2{t)) , (4.1^ 



= a2e-'-^y{t - r)z2{t - r) - d,Z2{t), 

dt 

Using simple algebraic manipulations, we get four feasible equilibria for the 
system fl4.1l) . namely, 

(a) trivial equilibrium Eq{^, 0, 0); 

(b) boundary equilibrium Ei (ai/61, 0, 0); 

(c) planner equilibrium i?2(a;,|/, 0) exists only when (HI) aiUi > bidi, 

(d) interior equilibrium E3{x*,y*,Z2) exists if (H2) aiaia2 > A. 
Where 

_ di _ oiOi — bidi ^ aia2 — Cid^e'^^'^ ^ d^e'^'^'^ 

! y ; 7 5 y ) 

«! CiOi Oia2 0.2 

^ _ aiaia2 — A 
^2 ~ 7 ■ 

OiC2a2 

Further, (H2) implies that 



A = cid^aie'^^^ + hidia2- 



1 J f aiaia2 — hidia2\ 

The characteristic equation for trivial equilibrium Eq{0, 0, 0) is given by 

{X - ai){X + di){X - ds) = 0. (4.2) 



The characteristic equation fl4.2p has one positive and two negative roots, 
hence equihbrium Eq is a unstable saddle point. 

Similarly, the characteristic equation for boundary equilibrium Ei is as 
follows: 

(A + ds){\ + ai) (\ + di- ^] = 0. (4.3) 



Clearly, all the eigenvalues are negative only when aiOi < bidi, which stabilize 
El, otherwise it is unstable. Again, the characteristic equation for planner 
equilibrium E2{x,y,0) becomes 

(A + 4 - a2ye~'^^^) (A^ + hxX + ciaixy) = 0. (4.4) 

Since, both roots of the quadratic equation A^ + bixX + ciaixy = have 
negative real parts, hence the equilibrium E2 is locally asymptotically stable 
if a2?/e"'^2T < fQj^ oil ^ > (1/^2) loga22//(i3 := Tcr- 

Finally, the characteristic equation for interior equilibrium -^3(0:*, y*, z^) is 
given as: 

A=^ + Ai(r)A' + A2{t)X + Asir) + {Bi{t)\^ + B2{t)\ + B3) e"^" = 0, (4.5) 
where 

Ai{t) = biX*{T) + da, A2(r) = bidsX*{T) + cia^x* {T)y* (r) , 

A^ir) = cidsaiX*{T)y*{r), = -a2y*{T)e-''^\ 

B2{r) = a2y*{T) {c2Z*2{r) - 6ix*(r)) e''^^^ 

B,{t) = a2X*{T)y*{T) {biC2z;{r) - Ciai/(r)) e'^^^ 

We write Ai, Bi in place of Aj(r), -Bj(T) for i = 1, 2, 3 in the rest of the analysis. 
The characteristic equation (14. 5 p can be rewritten as: 

P(A,r) + g(A,r)e-^" = 0, (4.6) 

where 

P(A, r) = A=^ + AiA' + A2A + A3, g(A, r) = BiX^ + B2X + B^. (4.7) 

When r = 0, the characteristic equation (14.51) becomes 

A^ + 6ix*A^ + {ciaix*y* + C2a2y*zl)\ + biC2a2X* y* zl = 0. (4.8) 

Since (A(0) + Ei(0))(A2(0) + ^2(0)) - (^3(0) + ^3(0)) = 6iCi«ix*V > 0, 
therefore, using Routh-Hurwitz criterion, all the solutions of the characteristic 
equation (14.81) have negative real parts. Thus the interior equilibrium E^ is 
locally asymptotically stable for r = if it exists. 

In the following, we investigate the existence of purely imaginary roots 
A = ici;(ci; > 0) of character istic equation (14.61) . We apply Beretta and Kuang 
(IBeretta and Kuangl . |2002| ) geometric criterion which gives the existence of 
purely imaginary roots of a characteristic equation with delay dependent co- 
efficients. 



Lemma 4.1. // (112) holds, then the following are satisfies: 

1. P(0,r) + g(0,r)^0, 

2. P{iuj, r) + Q{iuj, t) ^ for all u e R, 

3. limsup{|Q(A,r)/P(A,r)| : |A| ^ oo, i?eA > 0} < 1, 

4. F{uj, t) = \P{iuj, r) p — \Q{iuj, r) p for each r has at most a finite number 
of real zeros, 

5. each positive root uj{t) of F{u, t) = is continuous and differentiable in 
T whenever it exists. 

Proof 1. For r G [0,f), 

P(0, r) + g(0, r) = biC2a2X*y*z*e-'''^ ^ 0. 

2. P{iuj, r) + Q{iuj, r) = -{Ai + Bi)uj^ + A3 + B3 + i[-uj^ + {A2 + ^2)^] ^ 0. 

3. Since P(A,r) is a third degree polynomial in A and Q{\,t) second degree, 
hence, lim sup{|Q(A,r)/P(A,r)| : |A| ^ 00, PeA > 0} = < 1. 

4. Let F be defined as 

F{u,T) = \P{zu,T)\^-\Q{zco,r)\\ 

From 

\P{ico, r)|2 = a;6 ^ {Aj - 2A2) co^ + {Aj - + Al 

and 

\Q{zu, r)|2 = Bfu^ + _ 2P1P3) a;2 + bI 
we have 

P(a;, t) = uj^ + p{t)uj^ + g(r)a;^ + r(r), 

where 

_ (aia2 - Cid^e'^'^'^){aibia2 - (61 + 2ai)cici3e''2T) 

01^2 

/ \ 2 2 *2 *2 I 01, * * * —dor 2 *2 * — 2d9Tf2foi a;*+C92^) 

g(r) = C;^a]^x ?/ + 2oiC2a2a; y -^2^ — C20i2y z^e ^ -1-22-1^ 
r(r) = — Scidsai + (aitti — b\di)a2e~'^'^'^ ■ 
It is obvious that property (iv) is satisfied. 

5. Since P(cj, r) is continuous in lo and r and it is differentiable with respect 
to bj. Therefore, from Implicit Function Theorem each root of P(w, r) = is 
continuous and differentiable in r. 

Hence all the conditions of the Lemma are satisfied, which ensure the 
existence of purely imaginary roots for the characteristic equation fl4.5p . □ 



Now let \ = iu {u > 0) he a. root of fl4.5l) . Substituting it into fl4.5p and 
separating the real and imaginary parts, we get 

B2U)smojT + (i?3 — -Bio;^) cos car = Aiu^ — A3, 



{Biu'^ — Bs) sin cur + B2UJCOSUJT = co^ — A2UJ, (4.9) 
which gives 

BiLO^ - {A1B2 + A2B1 + B3) LO^ + {A2B3 + A3B2) CO 



smooT 



Bluj^ + {Bl - 2BiB3)uj^ + Bl 



{B2 - Aifii) + {AiB3 + A3B, - A2B2) u? - A3B3 ^ ^. 
^°^""= B!u^ + (Bl - 2B,B3)u- + Bl • ('-'^^ 

We can define the angle 6{t) G [0, 27r], V r > as the solution of fl4.10p : 
Siw^ _ [^^^^ ^ ^^^^ + ^3 + [^^^^ + ^^5^] ^ 



sin 6(t) 



Blu^ + [5| - 2BiB3]lu^ + Bl 



^ , , [B2 - AiBi] a;^ + [AiBs + Agfii - ^2^2] cu^ - A3B3 

COS6{t) = ^ ; 1^ K ; — K , 4.11 

where u = uj{t) and such 6{t) is uniquely well defined for all r, so that 
F{uj{t),t) = 0. Hence 

UJ^ +p{T)uj^ + q{T)uj'^ + r{T) = 0. (4.12) 

Again the polynomial function F can be written as 

F(uJ,T) = h{cU^T), 

where h is a cubic polynomial, defined by 

h{z, t) := + p{t)z'^ + q{T)z + r{T). 

Applying the Descartes' rule of signs for the number of positive roots of 
h{z,T) = 0, we get the following four cases: 
Case I: Let 

/ii = {r > I p{t) > g(r) > and r(r) > 0}, 

lu = {r > \ rij) > 0, at least one p(r) < g(r) < 0}, 



h = hi U / 



12- 



In the interval J12, h{z,T) = either has or 2 positive roots. When the 
polynomial h{z, r) has no positive zero in I12, then it also has no positive zero 



in the interval Ji. Thus in this case, purely imaginary root of fl4.5p never exists. 
Case II: Let 

/21 = {r > I p{t) > and r(r) < 0}, 

= {r > I p{t) < 0, g(r) < and r(r) < 0}, 

/23 = {r > I P(r) < 0, g(r) > and r(r) < 0}. 

In the region I23, either one or three positive zeros of h{z, r) exist. Suppose 
only one positive zero is feasible in I23, then in the interval I2 = I21 U I22 U I23, 
h{z,T) has only one positive zero. Therefore, iu* with oj* = u}{t*) > is a 
purely imaginary root of f l4.5p if and only if r* is a zero of the Sn, where 



Sri 



9(T) + 2mi r ■ . 

T ^ , , , r G I2, with n G Nq. 

uj{t) 



Now we will verify the following lemma: 

Lemma 4.2 (Beretta and Kuang (jBeretta and Kuang . 20021 )). Assume that 
oj{t) is a positive real root of defined for r & I, I O R+q, and at some 



T* e I, 

Snir* 



0, for some n ^ Nq 



(4.13) 



Then a pair of simple conjugate pure imaginary roots A4.(r*) = iu{T*), A_(r*) = 
— zcj(r*) of lli4-5\ ) exists at t = t* which crosses the imaginary axis from left to 
right if 6{t*) > and crosses the imaginary axis from right to left if 5{t*) < 0, 
where 



S{t*) = Sign 



f dRe{X) 

y dT 



X=iu){T*) 



(w(T*),r*) } Sign 



dSnir) 



dT 



(4.14) 



Since dF{uj,T)/duj\^^^^^,^ 
fore, from f l4.14p . we get 



[6^5 + 4p(^)^3 + 2g(r)a;]^^^(^.) > 0, there- 



S{t*) = Sign 



dRe{\) 



dr 



Sign 



X=iuj{T* ) 



dSJr) 



dr 



Here, we can easily observe that S'„(0) < and 5'„(r) > S'„+i(r) V r G I2, 
n & Nq. Thus, if Sq has no zero in I2, then the function Sn also have no 
zero in I2 and if the function Sn has positive zeros, denoted by for some 
T E I2, n E Nq, then without loss of generality, we may assume that 



dT 



^ with Sn{Ti) = 0. 



Applying similar logic as in ( iBeretta and Kuangl . 120021 ). stability switches oc- 
cur at the zeros of Sq{t), denoted by Tq, if (H2) holds. Let us assume that 



T* = min{r G I2 \ So{t) = 0} and r** = max{r G I2 \ So{t) = 0}. 



Using the Hopf bifurcation theorem for functional differential equation (iHale and Lunel 



I993I ). we can conclude the existence of Hopf bifurcation in the following the- 



orem: 

Theorem 4.3. Let (H2 ) hold. The local behavior of the system Ii2. 1\) at inte- 
rior equilibrium is described as: 

1. If the function Sq{t) has no positive zero in I2, then the interior equilib- 
rium E^{x* , y*,Z2) is locally asymptotically stable for all t > 0. 

2. If the function ^^(r) has at least positive zero in I2 for some n & Nq, then 

is locally asymptotically stable for r G [0,r*) U (r**,r] and unstable 
and a Hopf bifurcation occurs for t G (r*,r**), i.e., stability switches of 
stability-instability-stability occur. 

Case III: If in the interval /a = /12, h{z, r) = has two positive roots, denoted 
by Ui and ^2, we get following two sequences of functions on I^: 

S«(.) = r - ^iW±^ i S<^>(r) = r - „ ^ JV„. 

where Oi{t) and O^i^r) are the solutions of fl4.11l) when u = uji,uj2 respectively. 
Similarly, we can also obtain for t G I3 that Sn\o) < and S^\t) > 5'^^^(r) 
with n e No, k = 1,2. Thus if S^^\t) > S^^\t), then Si^\T) > sP{t) 
and hence stability switch depends on all real roots of 5*1^^ (r) = 0, otherwise 
the stability switches depend on roots of both Sn\T) = and Sn\T) = 0. 
Furthermore, we can also obtain the similar results stated in Theorem 14.31 
Case IV: If in the interval J4 = I23, h{z,T) has three positive zeros, we can 
obtain the parallel results as in case III. 



5 Direction and stability of Hopf bifurcation 



In the previous section, we obtained the conditions, under which system (14.11) 
undergoes Hopf bifurcation, taking maturation delay (r) as the critical pa- 
rameter. Using the normal form theory a nd ce nter manifold reduction as 



described in Hassard et al. (IHassard et al.l . Il98ll ). we can investigate the di 



rection of Hopf bifurcation and the properties of these bifurcating periodic 
solutions. Hence, we always assume that system (14.11) undergoes Hopf bifur- 
cations at the critical value r* of r and there exists a pair of pure imaginary 
roots, i.e., ±iu{T*) of the characteristic equation fl4.5l) . 



Using the Appendix \^ we can compute the following values: 



fr\\ ^ f oi |2 lfi'02p\ , 921 

cm = ^ ^9209n - 2\9ii\ " + 

_ Re{C,m 
Re{X'iT*)}' 

/32 = 2i?e{Ci(0)}, (5.1) 
^ /,n{Ci(0)} + /i2/^{A-(r-)} 

J 2 — ) 

UJ*T* 

which determine the behavior of bifurcating periodic solution in the center 
manifold at the critical value r*, i.e., fi2 determines the direction of the Hopf 
bifurcation: if /i2 > (/i2 < 0), then the Hopf bifurcation is supercritical 
(subcritical) and the bifurcating periodic solution exists for r > r* (r < r*); 
P2 determines the stability of the bifurcating periodic solution: the bifurcating 
periodic solution is stable (unstable) if /32 < {(32 > 0) and T2 determines the 
period of the bifurcating periodic solution: the period increases (decreases) if 
T2 > (T2 < 0). 



6 Numerical simulation 

To verify the previously established results, consider a three species crop-pest- 
natural enemy stage structured food chain model with the following parameter 
values: 

^ = x(t) (ai-x(t)-y(t)), 

"^^^^^ yit) {a,x{t) - rfi - 0.6z2{t)) , (6.1) 



dt 

1.3e-°-^"?/(t - T)z2{t - r) - 0.3z2{t). 



dZ2{t) _ ^ ^^_o.4r 



dt 

If we choose ai = 1, ai = 0.1 and di = 0.5, then we obtain that for the 
trivial equilibrium £"0(0, 0, 0), the characteristic equation has three eigenvalues, 
Ai = 1, A2 = —0.5 and A3 = 0.3. This verifies that trivial equilibrium of the 
system (16. ip is a unstable saddle point. Similarly, the characteristic equation 
for the boundary equilibrium i?i(l,0,0) has Ai = — 1,A2 = —0.3 and A3 = 
—0.2 eigenvalues and hence the boundary equilibrium is locally asymptotically 
stable (see Figure [T](a)). 

Now let ai = 2, di = 0.05 and ai = 1.2. Then the condition (HI) for the 
positivity of equilibrium E2 is satisfy and the characteristic equation for the 
equilibrium £'2(0.0416667, 1.95833, 0) of the system (16. ip has three eigenvalues, 
namely, Ai = -0.0208333 - 0.312222z, A2 = -0.0208333 + 0.312222i and 



Figure 1: Nonnegative equilibrium points are stable (a) boundary equilibrium 
El is stable for oi = 1, di = 0.5, ai = 0.1; (b) planner equilibrium E2 is 
stable for ai = 2, di = 0.05, ai = 1.2 and r = 5.4 > Tcr = 5.34608. 




Figure 2: The interior equilibrium E3 is locally asymptotically stable for r = 0. 

A3 = —0.3 + 2.54583e~°'^^ eigenvalues. Here, A3 is negative only when r > 
f = 5.34608, in this case, the planner equilibrium E2 is locally asymptotically 
stable for r > 5.34608 (see Figure [Uj^b)). We, mainly focus on the dynamics 
of the interior equilibrium. To study the local behavior of interior equilibrium 
E3{x*{t), y*{T), (t)) of the system (16. ip . we take the same parameters values 
for tti, di and ai. We obtain that the equilibrium is positive if r < ri = 
5.34608 ( see Figure ED . 

To apply the Descartes' rule of signs, we plot the coefficients p(r), g(r) and 
r(r) of the function h{z,T) with the maturation delay r (see Figure H]). We 
can easily check that in the intervals /i2,-?^23, the function h{z,T) = has 
and 1 positive roots respectively. Thus exactly one zero of h{z, r) is feasible in 
only the interval I2 = [0, 2.59955) (see Figure|l](d)). Further, all the conditions 
of the lemma 14.11 are satisfied in the interval l2- Now taking 5*0 on one axis 
and r on another axis, we obtain that the function So{t) = t — 6{t)/uj{t) 
has two zeros r* ~ 0.743 and r** ^ 1.568 in the interval I2, i.e., there are 
two critical values of the maturation delay of natural enemies at which the 
stability switching occurs (see Figure E]). 
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Figure 3: The equilibrium E3{x* , y* , Z2) is positive for r < r = 5.34608. 
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Figure 4: (a) p(r); (b) g(r) and (c) r(r) in the interval [0,5.34608); (d) curve 
of a;(r) on r e [0,2.59955). 




Figure 5: (a) Graphs of functions 5*0 and 5*1 for r G [0,2.59955); (b) curve of 
5*0 for the same value of r. 

Using the same set of parametric values, it is obtained that the interior 
equilibrium E^i^x* , y* , Z2) is locally asymptotically stable if r < r* = 0.743 
and a Hopf bifurcation occurs if r > 0.743, see Figures [6t^a) and[6]|^a). The 
equilibrium again becomes locally asymptotically when r > r** = 1.568, see 
Figure [71 Thus the interior equilibrium of the system (16. ip is locally asymptot- 
ically stable for r e [0, 0.743)U(1.568, 5.34608] and is unstable for r G (r*, t**). 
Hence the stability switches from stability-instability- stability occurs. This is 
the verification of the Theorem 14.31 

Furthermore, for the system (16.11) . it is clear that Re{ci{0))\r=T* = —3.9481 
< and Re{ci{0))\r=T** = 9.3706 > 0, according to the formula given in 
section |5l Therefore, Hopf bifurcation for the interior equilibrium at r* (resp. 
r**) is forward (resp. backward) and the bifurcating periodic solution on the 
center manifold are orbitally asymptotically stable (see Figure |8]). Finally, we 
observe with the following set of parameters: ai = 7, bi = 1, Ci = 1, C2 = 
0.5, di = 0.05, d2 = 0.6, 0^3 = 1.2, ai = 1.5, ^2 = 2, then system (12. ip has 
a complex dynamics of multiple bifurcation (i.e., chaos) for the maturation 
delay r = 1.5 (see Figure |9]). 



7 Conclusion 

In this paper, we have proposed a three species crop-pest-natural enemy food 
chain mathematical model with stage structure and maturation delay for the 
natural enemy. We have studied the local stability of four nonnegative equi- 
libria of the system (12. ip . It is found that the trivial equilibrium Eq{0,0,0) 
is always unstable; the boundary equilibrium _E'i(ai/6i, 0, 0) is locally asymp- 
totically stable if aitti < bidi. Further, the planner equilibrium E2{x,y,0) is 
locally asymptotically stable if aiai > bidi and r > (1/^2) log(a2|//(i3), oth- 
erwise, it is unstable. The interior equilibrium is locally asymptotically stable 
if r < r*, a Hopf bifurcation occurs in the interval r* < r < r** and if the 
maturation delay crossed the second critical r**, then the interior equilibrium 




Figure 6: The dynamics of the interior equihbrium for first critical value of 
maturation delay r: (a) is locally asymptotically stable for r = 0.742 < 
r* = 0.743; (b) occurrence of Hopf bifurcation at for r = 0.75 > r* = 0.743. 




Figure 7: The dynamics of the interior equilibrium for second critical value of 
maturation delay r: (a) occurrence of Hopf bifurcation at E^ for r = 1.56 < 
r** = 1.568; (b) E^ is locally asymptotically stable for r = 1.57 > r** = 1.568. 
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Figure 8: (a) The bifurcated periodic solution for the interior equilibrium 
is locally asymptotically stable at r* = 0.743; (b) The bifurcated periodic 
solution for the interior equilibrium is locally asymptotically stable at r** = 
1.568. 




Figure 9: The chaotic behavior of the system at the interior equihbrium for 
parameter set: Oi = 7, bi = 1, Ci = 1, C2 = 0.5, di = 0.05, ^2 = 0.6, ^3 = 
1.2, ai = 1.5, ^2 = 2 and r = 1.5. 



becomes again stable. Thus the maturation delay plays an important role in 
switching of stability from stability-instability-stability. Furthermore, using a 
numerical simulation, we obtained that the existence of bifurcation and also 
observed the chaotic behavior of the system for a particular range of the mat- 
uration delay. In particular, it is observed that the larger maturation delay 
may lead to extinction of the natural enemy, i.e., natural enemies may extinct 
due to its stage structure. This shows that the maturation delay of natural en- 
emy is the controlling parameter for the pest population. Thus stage structure 
has a great importance in the dynamics of the three species crop-pest-natural 
enemy food chain. 



A Appendix 

Let, xi = X — X* , X2 = y — y*, x^ = Z2 — z^, Xi(t) = Xi^rt), r = r* -\- fi and 
dropping the bars for simplification of notations, system (14. ip is transformed 
into an FDE in C = C{[-1, 0], R^) as 



x{t) = L^{xt) + F{n,xt), (A.l) 

where x{t) = {xi{t) , X2(t) , X3{t))'^ e and : C ^ R, F : C x R ^ R are 
given, respectively, by 

-bix* -cix* \ / 0i(O) 
= (r*+/i) I a,y* -C2y* M^) 
-ds J \ 03(0) 

° \ / '^i(-i) 

(r*+/i) I U-i) 

tta^e-'^^Cr^+M) ^^^*g-d2(r*+M) J y 03(-l) 

(A.2) 



and 

F{fi, 0) = (r* + /i) ai0i(O)02(O) - C2(^2(O)03(O) | , (A.3) 
V «2e-'^^(^*+'^V2(-l)03(-l) 

where = (0i, 02, ^s)^ ^ C*- By the Riesz representation theorem, there exists 
a function t]{6, fi) of bounded variation for 6* G [—1, 0] such that 

L^0 = J r]{9,fi)^{9) ioiOeC. (A.4) 
In fact, we can choose 

-bix* —cix* \ 
r/(^,/i) = (r*+/i) I air -c^y* 6(6) 

-da / 

\ 

-(r* + /i) I + 

a2;22e"'^'^^*+''^ a2l/*e-'^2(r*+M) y 

(A.5) 

where 6 is Dirac deha function. For G C^([— 1, 0], i?^), define 

A(/i)0= ; (A.6) 

I J.idv{s,fi)(f){s), 9 = 0, 

and 

r 0, [-1,0), 
i?(/x)0 = <^ (A.7) 

Then the system (1A.1|) is equivalent to 

it = v4(/i)xt + (A.8) 
where Xt{e) = x{t + 9) for ^ G [-1, 0]. For i) G C^([0, 1], (i?^)*), define 

f ^e(0,i], 
= <^ (A.9) 

and a bihnear inner product 

{^{s),m) = mm- f f i^{i-e)dr^{e)mdi, (a.io) 

where ?7(^) = 77(^,0). Then A{0) and A* are adjoint operators. By the discus- 
sion in section m we know that ±iu}*T* are eigenvalues of ^(0). Thus, they 



are also eigenvalues of A*. We need to compute the eigenvector of A{0) and 
A* corresponding to iuj*T* and —iu*T*, respectively. 

Suppose that q{6) = (1, a, /3)^e*^'^*'^* is the eigenvector of ^4(0) correspond- 
ing to iuj*T*, then A{0)q{9) = iu*T*q{9). It follows from the definition of ^4(0) 
and ri{9, /i) that 



iu* + bix* 
-aiy* 




cix* 
iu* C2y* 
-a2Z2e~'^'^'^* e-iuj*T* ioj* - a2y*e~ 



g(0) = 0. 



e—iujT 



Then, we can easily obtain 
g(0) = (l,a,/3f, 

where a = -^^^^ and /3 = • 
Similarly, let q*{9) = D{1, a*, /3*)e~*^'^*'^* be the eigenvector of A* correspond- 
ing to —iuj*T*, then similarly we can obtain 



a 



bix* — iu* 



OiiV 

By flAlOj) we get 

{q*{s),q{9)) - 



(3* 



aiy*{iui* — + a2y*e 



pe _ 

/ D{1, a*, l3*)e-^*^'^^-'Ur]{9){l, a, /3) V"*"*^^^ 
l + aa* + /3/3*- (1, / (j){9)dr]{9){l,a, 



D 



D[l + aa* + + r* {az^ + (5y*) a2(5*e 



* — f/2T* —iu^r" 



Then we choose 
D 



p—d2T* p—iu)*T* 



l + aa* + 1313* + r* {az^ + a2l3*e 

such that g(^)) = 1 and = 0. 

In the following, we use the ideas in Hassard et al. ( iHassard et al.l . Il98ll ) 
to compute the coordinates describing center manifold Cq at /i = 0. Define 



z{t) = (g*, xt), W{t, 9) = Xt{Q) - 2Re[z{t)q 
On the center manifold Cq, we have 



(A.ll) 



z^ 



W{t,9) = W{z{t),zit).0) = W2o{9)- + Wni9)zz + Wo2i9)- + - ■ ■ , (A.12) 



where z and z are local coordinates for Cq in C in the direction of q* and q* . 
Note that W is real if Xj is real. We deal only with the real solution. For 
solution xt G Co of flA.8|l . since /i = 0, we have 

z{t) = iuj*T*z + q*{0)F{0,W{z,z,0) + 2Re{zq{0)}) 
=^ iu*T* + q*{0)Fo{z, z) = iuj*r* + g{z, z), 



where 



2 —2 

g{z, z) = r (O)Fo(z, ^) = ^20 + + 9o2{e)^ + ■■■ . (A.13) 



From (lA.llD and (IA.12p . we have 



and 

qi9) = il,a,f3fe^'-'^\ 
Thus, we can easily obtain that 

xu{0) = Wi^iof-^ + Wil\0)z-z + IVo?(0)y + z + z + O {\{zrz)\') , 

2 —2 

X2t{0) = Wil\0)j + Wil\0)zz + Wt,l\0)j + az + az + 0{\{z,z)\') 
^3t{0) = W^l\0)j + Wil\0)zz + W^l\0)j + pz + J-z + 0{\{z, z)\') 

X2t{-l)=wg\-lf-^ + wil\-l)z-z + + 



3\ 



From the definition of F{fi,Xt), we have 



-hxUO) - CiXit{0)x2t{0) 

g{z,z) = T*D{l,a*,/3*) \ aiXit(0)x2t(0) - C2X2t(0)a;3t(0) 

a2e-''2"*X2t(-l)x3t(-l) 

T*D {z"^ [-bi - a{ci - aia*) - C2a(3a* + aa2/3P*e-'^^^* e'^''''^'] 

+2zz [-bi - (ci - aia*)Re{a} - C2Re{a(3} + a2(3*e-'^^^' Re{a(3}] 

+z^ [-bi - (ci - aia*)a - C2a*a(3 + a2(3*a(3e-'^^^' e^''^'^'] 



+ -z^-z 



2biWiQ\0) - AbiWil\0) 
(ci - aia*) (^1^20^(0) + 2aWil\0) + 2Wil\o) + 1^2? (0) 
C2a* (/9W^io^(0) + 2f3wS\0) + 2aWS\o) + 01^20^(0)) 



+a2l3*e 



.(3) 



(2), 



20 



-l) + 2Pe-'''''^'wil\-l) 



-d2T* ^-2iLJ*T*l 



Comparing the coefficients with (lA.lSp . we obtain 

g2o = 2T*D[—bi — (ci — — C2(yf3a* + aa2f3/3*e~ 

gn = 2T*D[-bi - (ci - aia*)Re{a} - C2Re{a/3} + a2/3*e"'^'^' Re{a/3}], 
go2 = 2r*Z)[— 6i — (ci — aidi*)a — C2a*af3 + a2f3*af3e~ 



-d2T* ^2iu* T* 1 



921 



T*D 



-6i (2WiQ\0) + 4Wil\0) 
-(ci - aia*) {aW^l\0) + 2aWil\o) + 2Wi:f{0) + wS\o) 
-C2a* (mioi^) + 2f3Wil\0) + 2aWif{0) + aWS\o) 



+ 2ae-''^*^'wi:f{-l) + ae''''^'wS\-l) 



In order to determine 5^21, we need to compute W2o{9) and Wu{9). From flA.81) 
and flA.lip . we have 



AW - 2i?e[r (O)Fog(^)], 9 E [-1, 0), 

AW -2Re[r{0)Foq{9)]+Fo, 9 = 0, 
AW + H{z,z,9), (A.14) 



W = Xt — zq — zq 

def 



where 



2 —2 

H{z, z, 9) = H2o{9f- + Hir{9)z-z + i/o2(^)y + 



(A. 15) 



Note that on center manifold Cq near the origin 

W = W,z + W,t 
Thus, we obtain 

{A - 2ioo*T*)W2o{e) = -H2o{e), AWnie) = -Hn{e). (A.16) 
Comparing the coefficient with (lA.lSp gives that 

H2o{0) = -g2oq{0) - go2m, H^M = -giiq{9) - ^iig(^). (A.17) 
From flA.16p . flA.lTp and the definition of A, we have 

W2o{e) = 2zu*T*W2o{e) + g2oq{0) + ^02g(^). 
Noting q{e) = g(0)e^'^*^*^ hence 

W2o{e) = ^g(0)e^-*-*^ + -%g(0)e"^^-*-'^ + E,e'^-'^*\ (A.18) 

where Ei = {e[^\ e[^\ E^^) G i?^ is a constant vector. Similarly, from ( ]A.16P 
and flA.17p . we obtain 

^ii(^) = -^me'-'^'' + ^me-'-'^'' + E2, (A.19) 

UJ*T id T 

where E2 = {E^\ e'^\ Ef'^) G is also a constant vector. 

In the following we shall find out Ei and E2. From the definition of A and 
( 1A.16I) . we can obtain 

j d7]{e)W2o{e) = 2iu*T*W2m-H2m. (A.20) 

and 

j dr]{9)Wu{0) = -Hii{0), (A.21) 

where 7]{6) = 77(0,6'). From flA.14p and flA.lSp . we have 

—bi — Cia 

H2o{Q) = -^72og(0) - g^m + 2r* I aa^ - Caa/? I , (A.22) 



and 



—hi — ciRe{a} 

HniO) = -gnqiO) - ^iig(O) + 2r* | aiRe{a} - C2Re{aP} | . (A.23) 

a2e-'^^^* Re{af3 



Substituting f[08D and f[09D into f[O0D and noticing that 
(^iw*r*J - e*"*"*''c/r/(^)^ g(0) = 0, 

and 

we obtain 



*'dr,{eyj g(0) = 0, 



—61 — cia 

aai — C2a/3 



which leads to 



2iuj* + hix* 
-aiy* 



CiX 

2iu* 



C2y 



El 



-a2^2e""''"*e-2*'^*"* a2?/*e-'^2"*e-2^^*"* 



—61 — cia 

aai — C2a/3 



a2a/9e-'^2r*g-2i<^v* 



It follows that 
^1 ^ 



2 ( -a^y* 




2iw* 





C2y* 



—aiz^e. ^ e P2 — ^2?/ e e 



X 



—61 — C\a 
aa\ — C2a(3 



a2a/3e 



-d2T* —2ioj* T* 



Similarly, substituting f[09l) and f[03l) into f[OB . we can get 



^2 = 2 







CiX 







C2y 



-a2zle"'^^^* ^2 - a2|/*e-^2r- 



—hi — ciRe{a} 
X I aiRe{a} — C2Re{a/3} 
a2e-'^^^* Re{a^} 
(A.24) 



Thus, we can determine W20 and Wu from flA.181) and flA.191) . 
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